home *** CD-ROM | disk | FTP | other *** search
/ Business Shareware / Business Shareware.iso / start / hobby / aa_51 / cmoon.c < prev    next >
Text File  |  1993-02-09  |  36KB  |  1,342 lines
  1. /* cmoon.c
  2.  *
  3.  * Expansions for the geocentric ecliptic longitude, latitude,
  4.  * and distance of the Moon referred to the mean equinox and
  5.  * ecliptic of date.
  6.  *
  7.  * This program extends the ELP2000-85 analytical Lunar theory to
  8.  * cover a 22,000-year interval centered at J2000.  Relatively few
  9.  * changes were required to achieve a precision of about 0.1 arc
  10.  * minute on this interval.  The Delaunay arguments l, l', D, F
  11.  * are expressed by polynomials of degree 7.  All existing
  12.  * oscillatory terms were retained but 31 new terms of higher
  13.  * degree in time were added.
  14.  *
  15.  * The fit is much better in the remote past and future if
  16.  * secular terms are included in the arguments of the oscillatory
  17.  * perturbations.  Such adjustments cannot easily be incorporated
  18.  * into the 1991 lunar tables.  In this program the traditional
  19.  * literal arguments are used instead, with mean elements adjusted
  20.  * for a best fit to the reference ephemeris.  Reference positions
  21.  * were taken from a numerically integrated ephemeris that very
  22.  * accurately duplicates JPL's DE200 model of the Moon and planets.
  23.  *
  24.  * This program omits many oscillatory terms from the analytical
  25.  * theory which, if they were included, would yield a much higher
  26.  * accuracy for modern dates.  Detailed statistics of the precision
  27.  * are given in the file cmp.doc.  Comparing at 400-day and some
  28.  * 4-day intervals over the period -8981 to +12920, the maximum
  29.  * errors noted were 7" longitude, 8" latitude, and 8 x 10^-8 au radius.
  30.  * The expressions used for precession in this comparision were
  31.  * those of Laskar.
  32.  * 
  33.  * The new coefficients were found by an unweighted least squares
  34.  * fit to the numerical ephemeris in the mentioned test interval.
  35.  * The approximation error increases rapidly outside this interval.
  36.  *
  37.  * Although this program approximates a very precise projection
  38.  * of DE200 into the past and future, the true physical accuracy
  39.  * of the ephemeris is limited by uncertainty in the coefficient
  40.  * of tidal acceleration of the Moon.  The magnitude of this
  41.  * uncertainty is itself uncertain, but it is believed to lie
  42.  * between 0 and 1 arcsec per century squared in the Moon's
  43.  * longitude.
  44.  *
  45.  *
  46.  *
  47.  * References:
  48.  *
  49.  *   P. Bretagnon and Francou, G., "Planetary theories in rectangular
  50.  *   and spherical variables. VSOP87 solutions," Astronomy and
  51.  *   Astrophysics 202, 309-315 (1988)
  52.  *
  53.  *   M. Chapront-Touze' and J. Chapront, "ELP2000-85: a semi-analytical
  54.  *   lunar ephemeris adequate for historical times," Astronomy and
  55.  *   Astrophysics 190, 342-352 (1988).
  56.  *
  57.  *   M. Chapront-Touze' and J. Chapront, _Lunar Tables and
  58.  *   Programs from 4000 B.C. to A.D. 8000_, Willmann-Bell (1991)
  59.  *
  60.  *   J. Laskar, "Secular terms of classical planetary theories
  61.  *   using the results of general theory," Astronomy and Astrophysics
  62.  *   157, 59070 (1986)
  63.  *
  64.  *   S. L. Moshier, "Comparison of a 7000-year lunar ephemeris
  65.  *   with analytical theory," Astronomy and Astrophysics 262,
  66.  *   613-616 (1992)
  67.  *
  68.  *
  69.  *
  70.  *
  71.  *
  72.  * The program has two subroutine entry points.
  73.  * Entry gmoon() returns the geometric position of the Moon
  74.  * relative to the Earth.  Its calling procedure is as follows:
  75.  *
  76.  * double JD;       input Julian Ephemeris Date
  77.  * double rect[3];  output equatorial x, y, z, in astronomical units
  78.  * double pol[3];   output ecliptic polar coordinatees in radians and au
  79.  *                  pol[0] longitude, pol[1] latitude, pol[2] radius
  80.  *  gmoon( JD, rect, pol );
  81.  *
  82.  * The second entry, domoon(), is intended to be called by the AA.ARC
  83.  * ephemeris calculator.  It calculates the Moon's apparent position.
  84.  * Corrections for nutation and light time are included, but the
  85.  * required routines for this are given as separate files in AA.ARC.
  86.  * Note that the difference between Ephemeris and Universal
  87.  * time is significant here, since the Moon moves in right
  88.  * ascension at the rate of about 2s/minute.
  89.  *
  90.  *
  91.  * - S. L. Moshier, August, 1991
  92.  * Last rev: July, 1992
  93.  */
  94.  
  95.  
  96. /* Define 1 to generate AA.ARC subroutines, else define 0 */
  97. #define AASUB 1
  98.  
  99.  
  100. /* Perturbation tables
  101.  */
  102. #define NLR 118
  103. static short LR[8*NLR] = {
  104. /*
  105.                Longitude    Radius
  106.  D  l' l  F    1"  .0001"  1km  .0001km */
  107.  
  108.  0, 0, 1, 0, 22639, 5858,-20905,-3550,
  109.  2, 0,-1, 0,  4586, 4383, -3699,-1109,
  110.  2, 0, 0, 0,  2369, 9139, -2955,-9676,
  111.  0, 0, 2, 0,   769,  257,  -569,-9251,
  112.  0, 1, 0, 0,  -666,-4171,    48, 8883,
  113.  0, 0, 0, 2,  -411,-5957,    -3,-1483,
  114.  2, 0,-2, 0,   211, 6556,   246, 1585,
  115.  2,-1,-1, 0,   205, 4358,  -152,-1377,
  116.  2, 0, 1, 0,   191, 9562,  -170,-7331,
  117.  2,-1, 0, 0,   164, 7285,  -204,-5860,
  118.  0, 1,-1, 0,  -147,-3213,  -129,-6201,
  119.  1, 0, 0, 0,  -124,-9881,   108, 7427,
  120.  0, 1, 1, 0,  -109,-3803,   104, 7552,
  121.  2, 0, 0,-2,    55, 1771,    10, 3211,
  122.  0, 0, 1, 2,   -45, -996,     0,    0,
  123.  0, 0, 1,-2,    39, 5333,    79, 6606,
  124.  4, 0,-1, 0,    38, 4298,   -34,-7825,
  125.  0, 0, 3, 0,    36, 1238,   -23,-2104,
  126.  4, 0,-2, 0,    30, 7726,   -21,-6363,
  127.  2, 1,-1, 0,   -28,-3971,    24, 2085,
  128.  2, 1, 0, 0,   -24,-3582,    30, 8238,
  129.  1, 0,-1, 0,   -18,-5847,    -8,-3791,
  130.  1, 1, 0, 0,    17, 9545,   -16,-6747,
  131.  2,-1, 1, 0,    14, 5303,   -12,-8314,
  132.  2, 0, 2, 0,    14, 3797,   -10,-4448,
  133.  4, 0, 0, 0,    13, 8991,   -11,-6500,
  134.  2, 0,-3, 0,    13, 1941,    14, 4027,
  135.  0, 1,-2, 0,    -9,-6791,    -7,  -27,
  136.  2, 0,-1, 2,    -9,-3659,     0, 7740,
  137.  2,-1,-2, 0,     8, 6055,    10,  562,
  138.  1, 0, 1, 0,    -8,-4531,     6, 3220,
  139.  2,-2, 0, 0,     8,  502,    -9,-8845,
  140.  0, 1, 2, 0,    -7,-6302,     5, 7509,
  141.  0, 2, 0, 0,    -7,-4475,     1,  657,
  142.  2,-2,-1, 0,     7, 3712,    -4,-9501,
  143.  2, 0, 1,-2,    -6,-3832,     4, 1311,
  144.  2, 0, 0, 2,    -5,-7416,     0,    0,
  145.  4,-1,-1, 0,     4, 3740,    -3,-9580,
  146.  0, 0, 2, 2,    -3,-9976,     0,    0,
  147.  3, 0,-1, 0,    -3,-2097,     3, 2582,
  148.  2, 1, 1, 0,    -2,-9145,     2, 6164,
  149.  4,-1,-2, 0,     2, 7319,    -1,-8970,
  150.  0, 2,-1, 0,    -2,-5679,    -2,-1171,
  151.  2, 2,-1, 0,    -2,-5212,     2, 3536,
  152.  2, 1,-2, 0,     2, 4889,     0, 1437,
  153.  2,-1, 0,-2,     2, 1461,     0, 6571,
  154.  4, 0, 1, 0,     1, 9777,    -1,-4226,
  155.  0, 0, 4, 0,     1, 9337,    -1,-1169,
  156.  4,-1, 0, 0,     1, 8708,    -1,-5714,
  157.  1, 0,-2, 0,    -1,-7530,    -1,-7385,
  158.  2, 1, 0,-2,    -1,-4372,     0,-1357,
  159.  0, 0, 2,-2,    -1,-3726,    -4,-4212,
  160.  1, 1, 1, 0,     1, 2618,     0,-9333,
  161.  3, 0,-2, 0,    -1,-2241,     0, 8624,
  162.  4, 0,-3, 0,     1, 1868,     0,-5142,
  163.  2,-1, 2, 0,     1, 1770,     0,-8488,
  164.  0, 2, 1, 0,    -1,-1617,     1, 1655,
  165.  1, 1,-1, 0,     1,  777,     0, 8512,
  166.  2, 0, 3, 0,     1,  595,     0,-6697,
  167.  2, 0, 1, 2,     0,-9902,     0,    0,
  168.  2, 0,-4, 0,     0, 9483,     0, 7785,
  169.  2,-2, 1, 0,     0, 7517,     0,-6575,
  170.  0, 1,-3, 0,     0,-6694,     0,-4224,
  171.  4, 1,-1, 0,     0,-6352,     0, 5788,
  172.  1, 0, 2, 0,     0,-5840,     0, 3785,
  173.  1, 0, 0,-2,     0,-5833,     0,-7956,
  174.  6, 0,-2, 0,     0, 5716,     0,-4225,
  175.  2, 0,-2,-2,     0,-5606,     0, 4726,
  176.  1,-1, 0, 0,     0,-5569,     0, 4976,
  177.  0, 1, 3, 0,     0,-5459,     0, 3551,
  178.  2, 0,-2, 2,     0,-5357,     0, 7740,
  179.  2, 0,-1,-2,     0, 1790,     8, 7516,
  180.  3, 0, 0, 0,     0, 4042,    -1,-4189,
  181.  2,-1,-3, 0,     0, 4784,     0, 4950,
  182.  2,-1, 3, 0,     0,  932,     0, -585,
  183.  2, 0, 2,-2,     0,-4538,     0, 2840,
  184.  2,-1,-1, 2,     0,-4262,     0,  373,
  185.  0, 0, 0, 4,     0, 4203,     0,    0,
  186.  0, 1, 0, 2,     0, 4134,     0,-1580,
  187.  6, 0,-1, 0,     0, 3945,     0,-2866,
  188.  2,-1, 0, 2,     0,-3821,     0,    0,
  189.  2,-1, 1,-2,     0,-3745,     0, 2094,
  190.  4, 1,-2, 0,     0,-3576,     0, 2370,
  191.  1, 1,-2, 0,     0, 3497,     0, 3323,
  192.  2,-3, 0, 0,     0, 3398,     0,-4107,
  193.  0, 0, 3, 2,     0,-3286,     0,    0,
  194.  4,-2,-1, 0,     0,-3087,     0,-2790,
  195.  0, 1,-1,-2,     0, 3015,     0,    0,
  196.  4, 0,-1,-2,     0, 3009,     0,-3218,
  197.  2,-2,-2, 0,     0, 2942,     0, 3430,
  198.  6, 0,-3, 0,     0, 2925,     0,-1832,
  199.  2, 1, 2, 0,     0,-2902,     0, 2125,
  200.  4, 1, 0, 0,     0,-2891,     0, 2445,
  201.  4,-1, 1, 0,     0, 2825,     0,-2029,
  202.  3, 1,-1, 0,     0, 2737,     0,-2126,
  203.  0, 1, 1, 2,     0, 2634,     0,    0,
  204.  1, 0, 0, 2,     0, 2543,     0,    0,
  205.  3, 0, 0,-2,     0,-2530,     0, 2010,
  206.  2, 2,-2, 0,     0,-2499,     0,-1089,
  207.  2,-3,-1, 0,     0, 2469,     0,-1481,
  208.  3,-1,-1, 0,     0,-2314,     0, 2556,
  209.  4, 0, 2, 0,     0, 2185,     0,-1392,
  210.  4, 0,-1, 2,     0,-2013,     0, 0,
  211.  0, 2,-2, 0,     0,-1931,     0, 0,
  212.  2, 2, 0, 0,     0,-1858,     0, 0,
  213.  2, 1,-3, 0,     0, 1762,     0, 0,
  214.  4, 0,-2, 2,     0,-1698,     0, 0,
  215.  4,-2,-2, 0,     0, 1578,     0,-1083,
  216.  4,-2, 0, 0,     0, 1522,     0,-1281,
  217.  3, 1, 0, 0,     0, 1499,     0,-1077,
  218.  1,-1,-1, 0,     0,-1364,     0, 1141,
  219.  1,-3, 0, 0,     0,-1281,     0, 0,
  220.  6, 0, 0, 0,     0, 1261,     0, -859,
  221.  2, 0, 2, 2,     0,-1239,     0, 0,
  222.  1,-1, 1, 0,     0,-1207,     0, 1100,
  223.  0, 0, 5, 0,     0, 1110,     0, -589,
  224.  0, 3, 0, 0,     0,-1013,     0,  213,
  225.  4,-1,-3, 0,     0,  998,     0, 0,
  226. };
  227.  
  228. #define NMB 56
  229. static short MB[6*NMB] = {
  230. /*
  231.                Latitude
  232.  D  l' l  F    1"  .0001" */
  233.  
  234.  0, 0, 0, 1,18461, 2387,
  235.  0, 0, 1, 1, 1010, 1671,
  236.  0, 0, 1,-1,  999, 6936,
  237.  2, 0, 0,-1,  623, 6524,
  238.  2, 0,-1, 1,  199, 4837,
  239.  2, 0,-1,-1,  166, 5741,
  240.  2, 0, 0, 1,  117, 2607,
  241.  0, 0, 2, 1,   61, 9120,
  242.  2, 0, 1,-1,   33, 3572,
  243.  0, 0, 2,-1,   31, 7597,
  244.  2,-1, 0,-1,   29, 5766,
  245.  2, 0,-2,-1,   15, 5663,
  246.  2, 0, 1, 1,   15, 1216,
  247.  2, 1, 0,-1,  -12, -941,
  248.  2,-1,-1, 1,    8, 8681,
  249.  2,-1, 0, 1,    7, 9586,
  250.  2,-1,-1,-1,    7, 4346,
  251.  0, 1,-1,-1,   -6,-7314,
  252.  4, 0,-1,-1,    6, 5796,
  253.  0, 1, 0, 1,   -6,-4601,
  254.  0, 0, 0, 3,   -6,-2965,
  255.  0, 1,-1, 1,   -5,-6324,
  256.  1, 0, 0, 1,   -5,-3684,
  257.  0, 1, 1, 1,   -5,-3113,
  258.  0, 1, 1,-1,   -5, -759,
  259.  0, 1, 0,-1,   -4,-8396,
  260.  1, 0, 0,-1,   -4,-8057,
  261.  0, 0, 3, 1,    3, 9841,
  262.  4, 0, 0,-1,    3, 6745,
  263.  4, 0,-1, 1,    2, 9985,
  264.  0, 0, 1,-3,    2, 7986,
  265.  4, 0,-2, 1,    2, 4139,
  266.  2, 0, 0,-3,    2, 1863,
  267.  2, 0, 2,-1,    2, 1462,
  268.  2,-1, 1,-1,    1, 7660,
  269.  2, 0,-2, 1,   -1,-6244,
  270.  0, 0, 3,-1,    1, 5813,
  271.  2, 0, 2, 1,    1, 5198,
  272.  2, 0,-3,-1,    1, 5156,
  273.  2, 1,-1, 1,   -1,-3178,
  274.  2, 1, 0, 1,   -1,-2643,
  275.  4, 0, 0, 1,    1, 1919,
  276.  2,-1, 1, 1,    1, 1346,
  277.  2,-2, 0,-1,    1,  859,
  278.  0, 0, 1, 3,   -1, -194,
  279.  2, 1, 1,-1,    0,-8227,
  280.  1, 1, 0,-1,    0, 8042,
  281.  1, 1, 0, 1,    0, 8026,
  282.  0, 1,-2,-1,    0,-7932,
  283.  2, 1,-1,-1,    0,-7910,
  284.  1, 0, 1, 1,    0,-6674,
  285.  2,-1,-2,-1,    0, 6502,
  286.  0, 1, 2, 1,    0,-6388,
  287.  4, 0,-2,-1,    0, 6337,
  288.  4,-1,-1,-1,    0, 5958,
  289.  1, 0, 1,-1,    0,-5889,
  290. };
  291.  
  292. #define NLRT 38
  293. static short LRT[8*NLRT] = {
  294. /*
  295. Multiply by T
  296.                Longitude    Radius
  297.  D  l' l  F   .1"  .00001" .1km  .00001km */
  298.  
  299.  0, 1, 0, 0,    16, 7680,    -1,-2302,
  300.  2,-1,-1, 0,    -5,-1642,     3, 8245,
  301.  2,-1, 0, 0,    -4,-1383,     5, 1395,
  302.  0, 1,-1, 0,     3, 7115,     3, 2654,
  303.  0, 1, 1, 0,     2, 7560,    -2,-6396,
  304.  2, 1,-1, 0,     0, 7118,     0,-6068,
  305.  2, 1, 0, 0,     0, 6128,     0,-7754,
  306.  1, 1, 0, 0,     0,-4516,     0, 4194,
  307.  2,-2, 0, 0,     0,-4048,     0, 4970,
  308.  0, 2, 0, 0,     0, 3747,     0, -540,
  309.  2,-2,-1, 0,     0,-3707,     0, 2490,
  310.  2,-1, 1, 0,     0,-3649,     0, 3222,
  311.  0, 1,-2, 0,     0, 2438,     0, 1760,
  312.  2,-1,-2, 0,     0,-2165,     0,-2530,
  313.  0, 1, 2, 0,     0, 1923,     0,-1450,
  314.  0, 2,-1, 0,     0, 1292,     0, 1070,
  315.  2, 2,-1, 0,     0, 1271,     0,-6070,
  316.  4,-1,-1, 0,     0,-1098,     0,  990,
  317.  2, 0, 0, 0,     0, 1073,     0,-1360,
  318.  2, 0,-1, 0,     0,  839,     0, -630,
  319.  2, 1, 1, 0,     0,  734,     0, -660,
  320.  4,-1,-2, 0,     0, -688,     0,  480,
  321.  2, 1,-2, 0,     0, -630,     0,    0,
  322.  0, 2, 1, 0,     0,  587,     0, -590,
  323.  2,-1, 0,-2,     0, -540,     0, -170,
  324.  4,-1, 0, 0,     0, -468,     0,  390,
  325.  2,-2, 1, 0,     0, -378,     0,  330,
  326.  2, 1, 0,-2,     0,  364,     0,    0,
  327.  1, 1, 1, 0,     0, -317,     0,  240,
  328.  2,-1, 2, 0,     0, -295,     0,  210,
  329.  1, 1,-1, 0,     0, -270,     0, -210,
  330.  2,-3, 0, 0,     0, -256,     0,  310,
  331.  2,-3,-1, 0,     0, -187,     0,  110,
  332.  0, 1,-3, 0,     0,  169,     0,  110,
  333.  4, 1,-1, 0,     0,  158,     0, -150,
  334.  4,-2,-1, 0,     0, -155,     0,  140,
  335.  0, 0, 1, 0,     0,  155,     0, -250,
  336.  2,-2,-2, 0,     0, -148,     0, -170,
  337. };
  338.  
  339. #define NBT 16
  340. static short BT[5*NBT] = {
  341. /*
  342. Multiply by T
  343.              Latitude
  344.  D  l' l  F  .00001"  */
  345.  
  346.  2,-1, 0,-1, -7430,
  347.  2, 1, 0,-1,  3043,
  348.  2,-1,-1, 1, -2229,
  349.  2,-1, 0, 1, -1999,
  350.  2,-1,-1,-1, -1869,
  351.  0, 1,-1,-1,  1696,
  352.  0, 1, 0, 1,  1623,
  353.  0, 1,-1, 1,  1418,
  354.  0, 1, 1, 1,  1339,
  355.  0, 1, 1,-1,  1278,
  356.  0, 1, 0,-1,  1217,
  357.  2,-2, 0,-1,  -547,
  358.  2,-1, 1,-1,  -443,
  359.  2, 1,-1, 1,   331,
  360.  2, 1, 0, 1,   317,
  361.  2, 0, 0,-1,   295,
  362. };
  363.  
  364. #define NLRT2 25
  365. static short LRT2[6*NLRT2] = {
  366. /*
  367. Multiply by T^2
  368.            Longitude    Radius
  369.  D  l' l  F  .00001" .00001km   */
  370.  
  371.  0, 1, 0, 0,  487,   -36,
  372.  2,-1,-1, 0, -150,   111,
  373.  2,-1, 0, 0, -120,   149,
  374.  0, 1,-1, 0,  108,    95,
  375.  0, 1, 1, 0,   80,   -77,
  376.  2, 1,-1, 0,   21,   -18,
  377.  2, 1, 0, 0,   20,   -23,
  378.  1, 1, 0, 0,  -13,    12,
  379.  2,-2, 0, 0,  -12,    14,
  380.  2,-1, 1, 0,  -11,     9,
  381.  2,-2,-1, 0,  -11,     7,
  382.  0, 2, 0, 0,   11,     0,
  383.  2,-1,-2, 0,   -6,    -7,
  384.  0, 1,-2, 0,    7,     5,
  385.  0, 1, 2, 0,    6,    -4,
  386.  2, 2,-1, 0,    5,    -3,
  387.  0, 2,-1, 0,    5,     3,
  388.  4,-1,-1, 0,   -3,     3,
  389.  2, 0, 0, 0,    3,    -4,
  390.  4,-1,-2, 0,   -2,     0,
  391.  2, 1,-2, 0,   -2,     0,
  392.  2,-1, 0,-2,   -2,     0,
  393.  2, 1, 1, 0,    2,    -2,
  394.  2, 0,-1, 0,    2,     0,
  395.  0, 2, 1, 0,    2,     0,
  396. };
  397.  
  398. #define NBT2 12
  399. static short BT2[5*NBT2] = {
  400. /*
  401. Multiply by T^2
  402.            Latitiude
  403.  D  l' l  F  .00001" */
  404.  
  405.  2,-1, 0,-1,  -22,
  406.  2, 1, 0,-1,    9,
  407.  2,-1, 0, 1,   -6,
  408.  2,-1,-1, 1,   -6,
  409.  2,-1,-1,-1,   -5,
  410.  0, 1, 0, 1,    5,
  411.  0, 1,-1,-1,    5,
  412.  0, 1, 1, 1,    4,
  413.  0, 1, 1,-1,    4,
  414.  0, 1, 0,-1,    4,
  415.  0, 1,-1, 1,    4,
  416.  2,-2, 0,-1,   -2,
  417. };
  418.  
  419. /* Rate of motion in R.A. and Dec., computed by this program.
  420.  */
  421. extern double dradt, ddecdt;
  422.  
  423. /* Obliquity of the ecliptic, computed by epsiln():
  424.  */
  425. extern double coseps, sineps;
  426.  
  427. /* Answers posted by angles():
  428.  */
  429. extern double pq, qe, ep;
  430.  
  431. /* Nutation, computed by nutlo():
  432.  */
  433. extern double nutl, nuto;
  434.  
  435. /* The following times are set up by update() and refer
  436.  * to the same instant.  The distinction between them
  437.  * is required by altaz().
  438.  */
  439. extern double TDT; /* Ephemeris time */
  440. extern double UT;  /* Universal time */
  441.  
  442. extern double ss[5][8]; /* see nutate.c */
  443. extern double cc[5][8];
  444.  
  445. #define DEBUG 0
  446. static double l = 0.0;        /* Moon's ecliptic longitude */
  447. static double B = 0.0;        /* Ecliptic latitude */
  448. static double p = 0.0;        /* Parallax */
  449. static double ra = 0.0;        /* Right Ascension */
  450. static double dec = 0.0;    /* Declination */
  451.  
  452. double moonpol[3];
  453. double moonpp[3];
  454.  
  455. /* Beware of atan2()!
  456.  */
  457. double floor(), sin(), cos(), asin(), mods3600();
  458.  
  459. /* equatorial radius of the earth, in au
  460.  */
  461. #define Rearth 4.263523e-5
  462. /* in kilometers */
  463. #define Kearth 6378.14
  464.  
  465. extern double STR, DTR, J2000;
  466.  
  467.  
  468.  
  469. #if AASUB
  470.  
  471. #include "kep.h"
  472. int moonll(), moon1(), moon2(), moon3(), moon4(), chewm(), sscc();
  473.  
  474. /* Calculate geometric position of the Moon and apply
  475.  * approximate corrections to find apparent position,
  476.  * phase of the Moon, etc. for AA.ARC.
  477.  */
  478. int domoon()
  479. {
  480. int i;
  481. double ra0, dec0, r;
  482. double x, y, z, lon0;
  483. double pp[3], qq[3];
  484. double acos();
  485.  
  486. /* Compute obliquity of the ecliptic, coseps, and sineps
  487.  */
  488. epsiln(TDT);
  489.  
  490. /* Run the orbit calculation twice, at two different times,
  491.  * in order to find the rate of change of R.A. and Dec.
  492.  */
  493.  
  494. /* Calculate for 0.001 day ago
  495.  */
  496. i = prtflg;
  497. prtflg = 0; /* disable display */
  498. moonll(TDT-0.001);
  499. lonlat( rearth, TDT, eapolar ); /* precess earth to date */
  500. prtflg = i; /* reenable display */
  501. ra0 = ra;
  502. dec0 = dec;
  503. lon0 = l;
  504.  
  505. /* Calculate for present instant.
  506.  */
  507. moonll(TDT);
  508.  
  509. /* The rates of change.  These are used by altaz() to
  510.  * correct the time of rising, transit, and setting.
  511.  */
  512. dradt = 1000.0*(ra-ra0);
  513. ddecdt = 1000.0*(dec-dec0);
  514.  
  515. /* Post the ecliptic longitude and latitude, in radians,
  516.  * and the radius in au.
  517.  */
  518. obpolar[0] = l;
  519. obpolar[1] = B;
  520. r = 1.0/sin(p); /* distance in earth-radii */
  521. ra0 = Rearth*r; /* factor is radius of earth in au */
  522. obpolar[2] = ra0;
  523.  
  524. /* Rate of change in longitude, degrees per day
  525.  * used for phase of the moon
  526.  */
  527. lon0 = 1000.0*RTD*(l - lon0);
  528.  
  529. /* convert to ecliptic rectangular coordinates */
  530.  
  531. z = ra0 * cos(B);
  532. x = z * cos(l);
  533. y = z * sin(l);
  534. z = ra0 * sin(B);
  535. /* convert to equatorial coordinates */
  536. pp[0] = x;
  537. pp[1] = y * coseps - z * sineps;
  538. pp[2] = y * sineps + z * coseps;
  539.  
  540. /* Find sun-moon-earth angles */
  541. precess( rearth, TDT, -1 );
  542. for( i=0; i<3; i++ )
  543.     qq[i] = rearth[i] + pp[i];
  544. angles( pp, qq, rearth );
  545.  
  546. /* Display answers
  547.  */
  548. if( prtflg )
  549.     {
  550.     printf( "Apparent geocentric longitude %.3lf deg", RTD*(l+nutl) );
  551.     printf( "   latitude %.3lf deg\n", RTD*B );
  552.     printf( "Distance %.5lf Earth-radii\n", r );
  553.     printf( "Horizontal parallax" );
  554.     dms( p );
  555.     printf( "Semidiameter" );
  556.     x = 0.272453 * p  +  0.0799/RTS; /* AA page L6 */
  557.     dms( x );
  558.  
  559.     x = RTD * acos(-ep);
  560.     printf( "\nElongation from sun %.2lf deg,", x );
  561.     x = 0.5 * (1.0 + pq);
  562.     printf( "  Illuminated fraction %.2lf\n", x );
  563.  
  564. /* Find phase of the Moon by comparing Moon's longitude
  565.  * with Earth's longitude.
  566.  *
  567.  * The number of days before or past indicated phase is
  568.  * estimated by assuming the true longitudes change linearly
  569.  * with time.  These rates are estimated for the date, but
  570.  * do not stay constant.  The error can exceed 0.15 day in 4 days.
  571.  */
  572. /* Apparent longitude of sun.
  573.  * Moon's longitude was already corrected for light time.
  574.  */
  575.     y = eapolar[0] - 20.496/(RTS*eapolar[2]);
  576.     x = obpolar[0] - y;
  577.     x = modtp( x ) * RTD;    /* difference in longitude */
  578.     i = x/90;        /* number of quarters */
  579.     x = (x - i*90.0);    /* phase angle mod 90 degrees */
  580.  
  581. /* days per degree of phase angle */
  582.     z = 1.0/(lon0 - (0.9856/eapolar[2]));
  583.  
  584.     if( x > 45.0 )
  585.         {
  586.         y = -(x - 90.0)*z;
  587.         if( y > 1.0 )
  588.             printf( "Phase %.1lf days before ", y );
  589.         else
  590.             printf( "Phase %.2lf days before ", y );
  591.         i = (i+1) & 3;
  592.         }
  593.     else
  594.         {
  595.         y = x*z;
  596.         if( y > 1.0 )
  597.             printf( "Phase %.1lf days past ", y );
  598.         else
  599.             printf( "Phase %.2lf days past ", y );
  600.         }
  601.  
  602.  
  603.     switch(i)
  604.         {
  605.         case 0: printf( "Full Moon\n" ); break;
  606.         case 1: printf( "Third Quarter\n" ); break;
  607.         case 2: printf( "New Moon\n" ); break;
  608.         case 3: printf( "First Quarter\n" ); break;
  609.         }
  610.     } /* if prtflg */
  611.  
  612. printf( "    Apparent:  R.A." );
  613. hms(ra);
  614. printf( "Declination" );
  615. dms(dec);
  616. printf( "\n" );
  617.  
  618. /* Compute and display topocentric position (altaz.c)
  619.  */
  620. pp[0] = ra;
  621. pp[1] = dec;
  622. pp[2] = r * Rearth;
  623. altaz( pp, UT );
  624. return(0);
  625. }
  626. #endif /* AASUB */
  627.  
  628.  
  629. /* Orbit calculation begins.
  630.  */
  631. static double f;
  632. static double g;
  633. static double LP;
  634. static double M;
  635. static double MP;
  636. static double D;
  637. static double NF;
  638. static double Ve;
  639. static double Ea;
  640. static double Ma;
  641. static double Ju;
  642. static double Sa;
  643. static double T;
  644. static double T2;
  645. static double T4;
  646. static double  cg, sg;
  647. static double l1, l2, l3, l4;
  648.  
  649. /* The following coefficients were calculated by a simultaneous least
  650.  * squares fit between the analytical theory and the continued DE200
  651.  * numerically integrated ephemeris from 9000 BC to 13000 AD.
  652.  * See references to the array z[] later on in the program.
  653.  * The 71 coefficients were estimated from 42,529 Lunar positions.
  654.  */
  655. static double z[] = {
  656. -1.225346551567e+001, /* F, t^2 */
  657. -1.096676093208e-003, /* F, t^3 */
  658. -2.165750777942e-006, /* F, t^4 */
  659. -2.790392351314e-009, /* F, t^5 */
  660.  4.189032191814e-011, /* F, t^6 */
  661.  4.474984866301e-013, /* F, t^7 */
  662.  3.239398410335e+001, /* l, t^2 */
  663.  5.185305877294e-002, /* l, t^3 */
  664. -2.536291235258e-004, /* l, t^4 */
  665. -2.506365935364e-008, /* l, t^5 */
  666.  3.452144225877e-011, /* l, t^6 */
  667. -1.755312760154e-012, /* l, t^7 */
  668. -5.870522364514e+000, /* D, t^2 */
  669.  6.493037519768e-003, /* D, t^3 */
  670. -3.702060118571e-005, /* D, t^4 */
  671.  2.560078201452e-009, /* D, t^5 */
  672.  2.555243317839e-011, /* D, t^6 */
  673. -3.207663637426e-013, /* D, t^7 */
  674. -4.776684245026e+000, /* L, t^2 */
  675.  6.580112707824e-003, /* L, t^3 */
  676. -6.073960534117e-005, /* L, t^4 */
  677. -1.024222633731e-008, /* L, t^5 */
  678.  2.235210987108e-010, /* L, t^6 */
  679.  7.200592540556e-014, /* L, t^7 */
  680. -8.552017636339e+001, /* t^2 cos(18V - 16E - l) */
  681. -2.055794304596e+002, /* t^2 sin(18V - 16E - l) */
  682. -1.097555241866e+000, /* t^3 cos(18V - 16E - l) */
  683.  5.219423171002e-001, /* t^3 sin(18V - 16E - l) */
  684.  2.088802640755e-003, /* t^4 cos(18V - 16E - l) */
  685.  4.616541527921e-003, /* t^4 sin(18V - 16E - l) */
  686.  4.794930645807e+000, /* t^2 cos(10V - 3E - l) */
  687. -4.595134364283e+001, /* t^2 sin(10V - 3E - l) */
  688. -6.659812174691e-002, /* t^3 cos(10V - 3E - l) */
  689. -2.570048828246e-001, /* t^3 sin(10V - 3E - l) */
  690.  6.229863046223e-004, /* t^4 cos(10V - 3E - l) */
  691.  5.504368344700e-003, /* t^4 sin(10V - 3E - l) */
  692. -3.084830597278e+000, /* t^2 cos(8V - 13E) */
  693. -1.000471012253e+001, /* t^2 sin(8V - 13E) */
  694.  6.590112074510e-002, /* t^3 cos(8V - 13E) */
  695. -3.212573348278e-003, /* t^3 sin(8V - 13E) */
  696.  5.409038312567e-004, /* t^4 cos(8V - 13E) */
  697.  1.293377988163e-003, /* t^4 sin(8V - 13E) */
  698.  2.311794636111e+001, /* t^2 cos(4E - 8M + 3J) */
  699. -3.157036220040e+000, /* t^2 sin(4E - 8M + 3J) */
  700. -3.019293162417e+000, /* t^2 cos(18V - 16E) */
  701. -9.211526858975e+000, /* t^2 sin(18V - 16E) */
  702. -4.993704215784e-002, /* t^3 cos(18V - 16E) */
  703.  2.991187525454e-002, /* t^3 sin(18V - 16E) */
  704. -3.827414182969e+000, /* t^2 cos(18V - 16E - 2l) */
  705. -9.891527703219e+000, /* t^2 sin(18V - 16E - 2l) */
  706. -5.322093802878e-002, /* t^3 cos(18V - 16E - 2l) */
  707.  3.164702647371e-002, /* t^3 sin(18V - 16E - 2l) */
  708.  7.713905234217e+000, /* t^2 cos(2J - 5S) */
  709. -6.077986950734e+000, /* t^3 sin(2J - 5S) */
  710. -1.278232501462e-001, /* t^2 cos(L - F) */
  711.  4.760967236383e-001, /* t^2 sin(L - F) */
  712. -6.759005756460e-001, /* t^3 sin(l') */
  713.  1.655727996357e-003, /* t^4 sin(l') */
  714.  1.646526117252e-001, /* t^3 sin(2D - l') */
  715. -4.167078100233e-004, /* t^4 sin(2D - l') */
  716.  2.067529538504e-001, /* t^3 sin(2D - l' - l) */
  717. -5.219127398748e-004, /* t^4 sin(2D - l' - l) */
  718. -1.526335222289e-001, /* t^3 sin(l' - l) */
  719. -1.120545131358e-001, /* t^3 sin(l' + l) */
  720.  4.619472391553e-002, /* t^3 sin(2D - 2l') */
  721.  4.863621236157e-004, /* t^4 sin(2D - 2l') */
  722. -4.280059182608e-002, /* t^3 sin(2l') */
  723. -4.328378207833e-004, /* t^4 sin(2l') */
  724. -8.371028286974e-003, /* t^3 sin(2D - l) */
  725.  4.089447328174e-002, /* t^3 sin(2D - 2l' - l) */
  726. -1.238363006354e-002, /* t^3 sin(2D + 2l' - l) */
  727. };
  728.  
  729.  
  730. #if AASUB
  731. /* Calculate apparent latitude, longitude, and horizontal parallax
  732.  * of the Moon at Julian date J.
  733.  */
  734. int moonll(J)
  735. double J;
  736. {
  737.  
  738. /* Note, the time scale is supposed to be TDB but the difference
  739.  * TDT - TDB is negligible compared to the truncation error
  740.  * of the trigonometric series.
  741.  */
  742. T = (J-J2000)/36525.0;
  743. T2 = T*T;
  744. T4 = T2*T2;
  745. /* Program is broken up to get it thru micro compiler */
  746. moon1();
  747. moon2();
  748. moon3();
  749. moon4(1);
  750. /* Correct for nutation at date TDT.
  751.  */
  752. nutate( TDT, moonpp );
  753.  
  754. ra = zatan2(moonpp[0],moonpp[1]);
  755. dec = asin(moonpp[2]);
  756. return(0);
  757. }
  758. #endif /* AASUB */
  759.  
  760. /* Calculate geometric coordinates of Moon
  761.  * without light time or nutation correction.
  762.  */
  763. int gmoon(J, rect, pol)
  764. double J;
  765. double rect[], pol[];
  766. {
  767. double r;
  768. int i;
  769.  
  770. epsiln(J);
  771. T = (J-J2000)/36525.0;
  772. T2 = T*T;
  773. T4 = T2*T2;
  774. moon1();
  775. moon2();
  776. moon3();
  777. moon4(0);
  778. r = moonpol[2];
  779. for( i=0; i<3; i++ )
  780.     {
  781.     rect[i] = moonpp[i] * r;
  782.     pol[i] = moonpol[i];
  783.     }
  784. return(0);
  785. }
  786.  
  787.  
  788.  
  789. /* Reduce arc seconds modulo 360 degrees
  790.  * answer in arc seconds
  791.  */
  792. double mods3600(x)
  793. double x;
  794. {
  795. double lx;
  796. double floor();
  797.  
  798. lx = x;
  799. lx = lx - 1296000.0 * floor( lx/1296000.0 );
  800. return( lx );
  801. }
  802.  
  803.  
  804. int moon1()
  805. {
  806. double a;
  807.  
  808. /* Mean anomaly of sun = l' (J. Laskar) */
  809. M =  mods3600(   129596581.038354 * T +  1287104.76154 );
  810. M += ((((((((
  811.   1.62e-20 * T
  812. - 1.0390e-17 ) * T
  813. - 3.83508e-15 ) * T
  814. + 4.237343e-13 ) * T
  815. + 8.8555011e-11 ) * T
  816. - 4.77258489e-8 ) * T
  817. - 1.1297037031e-5 ) * T
  818. + 1.4732069041e-4 ) * T
  819. - 0.552891801772 ) * T2;
  820.  
  821. /* Mean distance of moon from its ascending node = F */
  822. NF = mods3600( 1739527263.0983 * T + 335779.55755 );
  823. /* Mean anomaly of moon = l */
  824. MP = mods3600( 1717915923.4728 * T +  485868.28096 );
  825. /* Mean elongation of moon = D */
  826. D = mods3600( 1602961601.4603 * T + 1072260.73512 );
  827. /* Mean longitude of moon */
  828. LP = mods3600( 1732564372.83264 * T +  785939.95571 );
  829.  
  830. /* Higher degree secular terms found by least squares fit */
  831. NF += (((((z[5] *T+z[4] )*T + z[3] )*T + z[2] )*T + z[1] )*T + z[0] )*T2;
  832. MP += (((((z[11]*T+z[10])*T + z[9] )*T + z[8] )*T + z[7] )*T + z[6] )*T2;
  833. D  += (((((z[17]*T+z[16])*T + z[15])*T + z[14])*T + z[13])*T + z[12])*T2;
  834. LP += (((((z[23]*T+z[22])*T + z[21])*T + z[20])*T + z[19])*T + z[18])*T2;
  835.  
  836. /* sensitivity of mean elements
  837.  *    delta argument = scale factor times delta amplitude (arcsec)
  838.  * cos l  9.0019 = mean eccentricity
  839.  * cos 2D 43.6
  840.  * cos F  11.2 (latitude term)
  841.  */
  842.  
  843. /* Mean longitudes of planets (Laskar, Bretagnon) */
  844.  
  845. Ve = mods3600( 210664136.4335482 * T + 655127.283046 );
  846. Ve += ((((((((
  847.   -9.36e-023 * T
  848.  - 1.95e-20 ) * T
  849.  + 6.097e-18 ) * T
  850.  + 4.43201e-15 ) * T
  851.  + 2.509418e-13 ) * T
  852.  - 3.0622898e-10 ) * T
  853.  - 2.26602516e-9 ) * T
  854.  - 1.4244812531e-5 ) * T
  855.  + 0.005871373088 ) * T2;
  856.  
  857. Ea = mods3600( 129597742.26669231  * T +  361679.214649 );
  858. Ea += (((((((( -1.16e-22 * T
  859.  + 2.976e-19 ) * T
  860.  + 2.8460e-17 ) * T
  861.  - 1.08402e-14 ) * T
  862.  - 1.226182e-12 ) * T
  863.  + 1.7228268e-10 ) * T
  864.  + 1.515912254e-7 ) * T
  865.  + 8.863982531e-6 ) * T
  866.  - 2.0199859001e-2 ) * T2;
  867.  
  868. Ma = mods3600(  68905077.59284 * T + 1279559.78866 );
  869. Ma += (-1.043e-5*T + 9.38012e-3)*T2;
  870.  
  871. Ju = mods3600( 10925660.428608 * T +  123665.342120 );
  872. Ju += (1.543273e-5*T - 3.06037836351e-1)*T2;
  873.  
  874. Sa = mods3600( 4399609.65932 * T + 180278.89694 );
  875. Sa += (( 4.475946e-8*T - 6.874806E-5 ) * T + 7.56161437443E-1)*T2;
  876.  
  877. sscc( 0, STR*D, 6 );
  878. sscc( 1, STR*M,  4 );
  879. sscc( 2, STR*MP, 4 );
  880. sscc( 3, STR*NF, 4 );
  881.  
  882. moonpol[0] = 0.0;
  883. moonpol[1] = 0.0;
  884. moonpol[2] = 0.0;
  885.  
  886. /* terms in T^2, scale 1.0 = 10^-5" */
  887. chewm( LRT2, NLRT2, 4, 2, moonpol );
  888. chewm( BT2, NBT2, 4, 4, moonpol );
  889.  
  890. f = 18 * Ve - 16 * Ea;
  891.  
  892. g = STR*(f - MP );  /* 18V - 16E - l */
  893. cg = cos(g);
  894. sg = sin(g);
  895. l = 6.367278 * cg + 12.747036 * sg;  /* t^0 */
  896. l1 = 23123.70 * cg - 10570.02 * sg;  /* t^1 */
  897. l2 = z[24] * cg + z[25] * sg;        /* t^2 */
  898. l3 = z[26] * cg + z[27] * sg;        /* t^3 */
  899. l4 = z[28] * cg + z[29] * sg;        /* t^4 */
  900. moonpol[2] += 5.01 * cg + 2.72 * sg;
  901.  
  902. g = STR * (10.*Ve - 3.*Ea - MP);
  903. cg = cos(g);
  904. sg = sin(g);
  905. l += -0.253102 * cg + 0.503359 * sg;
  906. l1 += 1258.46 * cg + 707.29 * sg;
  907. l2 += z[30] * cg + z[31] * sg;
  908. l3 += z[32] * cg + z[33] * sg;
  909. l4 += z[34] * cg + z[35] * sg;
  910.  
  911. g = STR*(8.*Ve - 13.*Ea);
  912. cg = cos(g);
  913. sg = sin(g);
  914. l += -0.187231 * cg - 0.127481 * sg;
  915. l1 += -319.87 * cg - 18.34 * sg;
  916. l2 += z[36] * cg + z[37] * sg;
  917. l3 += z[38] * cg + z[39] * sg;
  918. l4 += z[40] * cg + z[41] * sg;
  919.  
  920. a = 4.0*Ea - 8.0*Ma + 3.0*Ju;
  921. g = STR * a;
  922. cg = cos(g);
  923. sg = sin(g);
  924. l += -0.866287 * cg + 0.248192 * sg;
  925. l1 += 41.87 * cg + 1053.97 * sg;
  926. l2 += z[42] * cg + z[43] * sg;
  927.  
  928. g = STR*(a - MP);
  929. cg = cos(g);
  930. sg = sin(g);
  931. l += -0.165009 * cg + 0.044176 * sg;
  932. l1 += 4.67 * cg + 201.55 * sg;
  933.  
  934.  
  935. g = STR*f;  /* 18V - 16E */
  936. cg = cos(g);
  937. sg = sin(g);
  938. l += 0.330401 * cg + 0.661362 * sg;
  939. l1 += 1202.67 * cg - 555.59 * sg;
  940. l2 += z[44] * cg + z[45] * sg;
  941. l3 += z[46] * cg + z[47] * sg;
  942.  
  943. g = STR*(f - 2.0*MP );  /* 18V - 16E - 2l */
  944. cg = cos(g);
  945. sg = sin(g);
  946. l += 0.352185 * cg + 0.705041 * sg;
  947. l1 += 1283.59 * cg - 586.43 * sg;
  948. l2 += z[48] * cg + z[49] * sg;
  949. l3 += z[50] * cg + z[51] * sg;
  950.  
  951. g = STR * (2.0*Ju - 5.0*Sa);
  952. cg = cos(g);
  953. sg = sin(g);
  954. l += -0.034700 * cg + 0.160041 * sg;
  955. l2 += z[52] * cg + z[53] * sg;
  956.  
  957. g = STR * (LP - NF);
  958. cg = cos(g);
  959. sg = sin(g);
  960. l += 0.000116 * cg + 7.063040 * sg;
  961. l1 +=  298.8 * sg;
  962. l2 += z[54] * cg + z[55] * sg;
  963.  
  964.  
  965. /* T^3 terms */
  966. sg = sin( STR * M );
  967. l3 +=  z[56] * sg;
  968. l4 +=  z[57] * sg;
  969.  
  970. g = STR * (2.0*D - M);
  971. sg = sin(g);
  972. cg = cos(g);
  973. l3 +=  z[58] * sg;
  974. l4 +=  z[59] * sg;
  975. moonpol[2] +=  -0.2655 * cg * T;
  976.  
  977. g = g - STR * MP;
  978. sg = sin(g);
  979. l3 +=  z[60] * sg;
  980. l4 +=  z[61] * sg;
  981.  
  982. g = STR * (M - MP);
  983. l3 +=  z[62] * sin( g );
  984. moonpol[2] +=  -0.1568 * cos( g ) * T;
  985.  
  986. g = STR * (M + MP);
  987. l3 +=  z[63] * sin( g );
  988. moonpol[2] +=  0.1309 * cos( g ) * T;
  989.  
  990. g = STR * 2.0 * (D - M);
  991. sg = sin(g);
  992. l3 +=  z[64] * sg;
  993. l4 +=  z[65] * sg;
  994.  
  995. g = STR * 2.0 * M;
  996. sg = sin(g);
  997. l3 +=  z[66] * sg;
  998. l4 +=  z[67] * sg;
  999.  
  1000. g = STR * (2.0*D - MP);
  1001. sg = sin(g);
  1002. l3 +=  z[68] * sg;
  1003.  
  1004. g = STR * (2.0*(D - M) - MP);
  1005. sg = sin(g);
  1006. l3 +=  z[69] * sg;
  1007.  
  1008. g = STR * (2.0*(D + M) - MP);
  1009. sg = sin(g);
  1010. cg = cos(g);
  1011. l3 +=  z[70] * sg;
  1012. moonpol[2] +=   0.5568 * cg * T;
  1013.  
  1014. l2 += moonpol[0];
  1015.  
  1016. g = STR*(2.0*D - M - MP);
  1017. moonpol[2] +=  -0.1910 * cos( g ) * T;
  1018.  
  1019.  
  1020. moonpol[1] *= T;
  1021. moonpol[2] *= T;
  1022.  
  1023. /* terms in T */
  1024. moonpol[0] = 0.0;
  1025. chewm( BT, NBT, 4, 4, moonpol );
  1026. chewm( LRT, NLRT, 4, 1, moonpol );
  1027. g = STR*(f - MP - NF - 2355767.6); /* 18V - 16E - l - F */
  1028. moonpol[1] +=  -1127. * sin(g);
  1029. g = STR*(f - MP + NF - 235353.6); /* 18V - 16E - l + F */
  1030. moonpol[1] +=  -1123. * sin(g);
  1031. g = STR*(Ea + D + 51987.6);
  1032. moonpol[1] +=  1303. * sin(g);
  1033. g = STR*LP;
  1034. moonpol[1] +=  342. * sin(g);
  1035.  
  1036.  
  1037. g = STR*(2.*Ve - 3.*Ea);
  1038. cg = cos(g);
  1039. sg = sin(g);
  1040. l +=  -0.343550 * cg - 0.000276 * sg;
  1041. l1 +=  105.90 * cg + 336.53 * sg;
  1042.  
  1043. g = STR*(f - 2.*D); /* 18V - 16E - 2D */
  1044. cg = cos(g);
  1045. sg = sin(g);
  1046. l += 0.074668 * cg + 0.149501 * sg;
  1047. l1 += 271.77 * cg - 124.20 * sg;
  1048.  
  1049. g = STR*(f - 2.*D - MP);
  1050. cg = cos(g);
  1051. sg = sin(g);
  1052. l += 0.073444 * cg + 0.147094 * sg;
  1053. l1 += 265.24 * cg - 121.16 * sg;
  1054.  
  1055. g = STR*(f + 2.*D - MP);
  1056. cg = cos(g);
  1057. sg = sin(g);
  1058. l += 0.072844 * cg + 0.145829 * sg;
  1059. l1 += 265.18 * cg - 121.29 * sg;
  1060.  
  1061. g = STR*(f + 2.*(D - MP));
  1062. cg = cos(g);
  1063. sg = sin(g);
  1064. l += 0.070201 * cg + 0.140542 * sg;
  1065. l1 += 255.36 * cg - 116.79 * sg;
  1066.  
  1067. g = STR*(Ea + D - NF);
  1068. cg = cos(g);
  1069. sg = sin(g);
  1070. l += 0.288209 * cg - 0.025901 * sg;
  1071. l1 += -63.51 * cg - 240.14 * sg;
  1072.  
  1073. g = STR*(2.*Ea - 3.*Ju + 2.*D - MP);
  1074. cg = cos(g);
  1075. sg = sin(g);
  1076. l += 0.077865 * cg + 0.438460 * sg;
  1077. l1 += 210.57 * cg + 124.84 * sg;
  1078.  
  1079. g = STR*(Ea - 2.*Ma);
  1080. cg = cos(g);
  1081. sg = sin(g);
  1082. l += -0.216579 * cg + 0.241702 * sg;
  1083. l1 += 197.67 * cg + 125.23 * sg;
  1084.  
  1085. g = STR*(a + MP);
  1086. cg = cos(g);
  1087. sg = sin(g);
  1088. l += -0.165009 * cg + 0.044176 * sg;
  1089. l1 += 4.67 * cg + 201.55 * sg;
  1090.  
  1091. g = STR*(a + 2.*D - MP);
  1092. cg = cos(g);
  1093. sg = sin(g);
  1094. l += -0.133533 * cg + 0.041116 * sg;
  1095. l1 +=  6.95 * cg + 187.07 * sg;
  1096.  
  1097. g = STR*(a - 2.*D + MP);
  1098. cg = cos(g);
  1099. sg = sin(g);
  1100. l += -0.133430 * cg + 0.041079 * sg;
  1101. l1 +=  6.28 * cg + 169.08 * sg;
  1102.  
  1103. g = STR*(3.*Ve - 4.*Ea);
  1104. cg = cos(g);
  1105. sg = sin(g);
  1106. l += -0.175074 * cg + 0.003035 * sg;
  1107. l1 +=  49.17 * cg + 150.57 * sg;
  1108.  
  1109. g = STR*(2.*(Ea + D - MP) - 3.*Ju + 213534.);
  1110. l1 +=  158.4 * sin(g);
  1111. l1 += moonpol[0];
  1112.  
  1113. a = 0.1 * T; /* set amplitude scale of 1.0 = 10^-4 arcsec */
  1114. moonpol[1] *= a;
  1115. moonpol[2] *= a;
  1116. return(0);
  1117. }
  1118.  
  1119.  
  1120.  
  1121. int moon2()
  1122. {
  1123. /* terms in T^0 */
  1124. g = STR*(2*(Ea-Ju+D)-MP+648431.172);
  1125. l += 1.14307 * sin(g);
  1126. g = STR*(Ve-Ea+648035.568);
  1127. l += 0.82155 * sin(g);
  1128. g = STR*(3*(Ve-Ea)+2*D-MP+647933.184);
  1129. l += 0.64371 * sin(g);
  1130. g = STR*(Ea-Ju+4424.04);
  1131. l += 0.63880 * sin(g);
  1132. g = STR*(LP + MP - NF + 4.68);
  1133. l += 0.49331 * sin(g);
  1134. g = STR*(LP - MP - NF + 4.68);
  1135. l += 0.4914 * sin(g);
  1136. g = STR*(LP+NF+2.52);
  1137. l += 0.36061 * sin(g);
  1138. g = STR*(2.*Ve - 2.*Ea + 736.2);
  1139. l += 0.30154 * sin(g);
  1140. g = STR*(2.*Ea - 3.*Ju + 2.*D - 2.*MP + 36138.2);
  1141. l += 0.28282 * sin(g);
  1142. g = STR*(2.*Ea - 2.*Ju + 2.*D - 2.*MP + 311.0);
  1143. l += 0.24516 * sin(g);
  1144. g = STR*(Ea - Ju - 2.*D + MP + 6275.88);
  1145. l += 0.21117 * sin(g);
  1146. g = STR*(2.*(Ea - Ma) - 846.36);
  1147. l += 0.19444 * sin(g);
  1148. g = STR*(2.*(Ea - Ju) + 1569.96);
  1149. l -= 0.18457 * sin(g);
  1150. g = STR*(2.*(Ea - Ju) - MP - 55.8);
  1151. l += 0.18256 * sin(g);
  1152. g = STR*(Ea - Ju - 2.*D + 6490.08);
  1153. l += 0.16499 * sin(g);
  1154. g = STR*(Ea - 2.*Ju - 212378.4);
  1155. l += 0.16427 * sin(g);
  1156. g = STR*(2.*(Ve - Ea - D) + MP + 1122.48);
  1157. l += 0.16088 * sin(g);
  1158. g = STR*(Ve - Ea - MP + 32.04);
  1159. l -= 0.15350 * sin(g);
  1160. g = STR*(Ea - Ju - MP + 4488.88);
  1161. l += 0.14346 * sin(g);
  1162. g = STR*(2.*(Ve - Ea + D) - MP - 8.64);
  1163. l += 0.13594 * sin(g);
  1164. g = STR*(2.*(Ve - Ea - D) + 1319.76);
  1165. l += 0.13432 * sin(g);
  1166. g = STR*(Ve - Ea - 2.*D + MP - 56.16);
  1167. l -= 0.13122 * sin(g);
  1168. g = STR*(Ve - Ea + MP + 54.36);
  1169. l -= 0.12722 * sin(g);
  1170. g = STR*(3.*(Ve - Ea) - MP + 433.8);
  1171. l += 0.12539 * sin(g);
  1172. g = STR*(Ea - Ju + MP + 4002.12);
  1173. l += 0.10994 * sin(g);
  1174. g = STR*(20.*Ve - 21.*Ea - 2.*D + MP - 317511.72);
  1175. l += 0.10652 * sin(g);
  1176. g = STR*(26.*Ve - 29.*Ea - MP + 270002.52);
  1177. l += 0.10490 * sin(g);
  1178. g = STR*(3.*Ve - 4.*Ea + D - MP - 322765.56);
  1179. l += 0.10386 * sin(g);
  1180.  
  1181.  
  1182. g = STR*(LP+648002.556);
  1183. B =  8.04508 * sin(g);
  1184. g = STR*(Ea+D+996048.252);
  1185. B += 1.51021 * sin(g);
  1186. g = STR*(f - MP + NF + 95554.332);
  1187. B += 0.63037 * sin(g);
  1188. g = STR*(f - MP - NF + 95553.792);
  1189. B += 0.63014 * sin(g);
  1190. g = STR*(LP - MP + 2.9);
  1191. B +=  0.45587 * sin(g);
  1192. g = STR*(LP + MP + 2.5);
  1193. B +=  -0.41573 * sin(g);
  1194. g = STR*(LP - 2.0*NF + 3.2);
  1195. B +=  0.32623 * sin(g);
  1196. g = STR*(LP - 2.0*D + 2.5);
  1197. B +=  0.29855 * sin(g);
  1198. return(0);
  1199. }
  1200.  
  1201.  
  1202.  
  1203. int moon3()
  1204. {
  1205. /* terms in T^0 */
  1206. moonpol[0] = 0.0;
  1207. chewm( LR, NLR, 4, 1, moonpol );
  1208. chewm( MB, NMB, 4, 3, moonpol );
  1209. l += (((l4 * T + l3) * T + l2) * T + l1) * T * 1.0e-5;
  1210. moonpol[0] = LP + l + 1.0e-4 * moonpol[0];
  1211. moonpol[1] = 1.0e-4 * moonpol[1] + B;
  1212. moonpol[2] = 1.0e-4 * moonpol[2] + 385000.52899; /* kilometers */
  1213. return(0);
  1214. }
  1215.  
  1216.  
  1217.  
  1218. /* Compute final ecliptic polar coordinates
  1219.  * and convert to equatorial rectangular coordinates
  1220.  */
  1221. int moon4(ltflag)
  1222. int ltflag;
  1223. {
  1224. double cosB, sinB, cosL, sinL, sp;
  1225.  
  1226. sp = Kearth/moonpol[2];
  1227. p = asin( sp );
  1228. moonpol[2] /= 1.4959787066e8; /* Kilometers per au */
  1229.  
  1230. l = STR * mods3600( moonpol[0] );
  1231.  
  1232. /* Light time correction to longitude,
  1233.  * about 0.7".
  1234.  */
  1235. if( ltflag )
  1236.     l -= DTR * 0.0118 * sp;
  1237. moonpol[0] = l;
  1238.  
  1239. B = STR * moonpol[1];
  1240. moonpol[1] = B;
  1241.  
  1242. /* convert to equatorial system of date */
  1243. cosB = cos(B);
  1244. sinB = sin(B);
  1245. cosL = cos(l);
  1246. sinL = sin(l);
  1247.  
  1248. moonpp[0] = cosB*cosL;
  1249. moonpp[1] = coseps*cosB*sinL - sineps*sinB;
  1250. moonpp[2] = sineps*cosB*sinL + coseps*sinB;
  1251. return(0);
  1252. }
  1253.  
  1254.  
  1255.  
  1256.  
  1257. /* Program to step through the perturbation table
  1258.  */
  1259. int chewm( p, nlines, nangles, typflg, ans )
  1260. register short *p;
  1261. int nlines, nangles;
  1262. int typflg;
  1263. double ans[];
  1264. {
  1265. int i, j, k, k1, m;
  1266. register double cu, su, cv, sv, f;
  1267.  
  1268.  
  1269. for( i=0; i<nlines; i++ )
  1270.     {
  1271.     k1 = 0;
  1272.     sv = 0.0;
  1273.     cv = 0.0;
  1274.     for( m=0; m<nangles; m++ )
  1275.         {
  1276.         j = *p++; /* multiple angle factor */
  1277.         if( j )
  1278.             {
  1279.             k = j;
  1280.             if( j < 0 )
  1281.                 k = -k; /* make angle factor > 0 */
  1282. /* sin, cos (k*angle) from lookup table */
  1283.             su = ss[m][k-1];
  1284.             cu = cc[m][k-1];
  1285.             if( j < 0 )
  1286.                 su = -su; /* negative angle factor */
  1287.             if( k1 == 0 )
  1288.                 {
  1289. /* Set sin, cos of first angle. */
  1290.                 sv = su;
  1291.                 cv = cu;
  1292.                 k1 = 1;
  1293.                 }
  1294.             else
  1295.                 {
  1296. /* Combine angles by trigonometry. */
  1297.                 f =  su*cv + cu*sv;
  1298.                 cv = cu*cv - su*sv;
  1299.                 sv = f;
  1300.                 }
  1301.             }
  1302.         }
  1303. /* Accumulate
  1304.  */
  1305.     switch( typflg )
  1306.         {
  1307. /* large longitude and radius */
  1308.         case 1:
  1309.             j = *p++;
  1310.             k = *p++;
  1311.             ans[0] += (10000.0 * j  + k) * sv;
  1312.             j = *p++;
  1313.             k = *p++;
  1314.             if( k )
  1315.                 {
  1316.                 ans[2] += (10000.0 * j  + k) * cv;
  1317.                 }
  1318.             break;
  1319. /* longitude and radius */
  1320.         case 2:
  1321.             j = *p++;
  1322.             k = *p++;
  1323.             ans[0] += j * sv;
  1324.             ans[2] += k * cv;
  1325.             break;
  1326. /* large latitude */
  1327.         case 3:
  1328.             j = *p++;
  1329.             k = *p++;
  1330.             ans[1] += ( 10000.0*j + k)*sv;
  1331.             break;
  1332. /* latitude */
  1333.         case 4:
  1334.             j = *p++;
  1335.             ans[1] += j * sv;
  1336.             break;
  1337.         }
  1338.     }
  1339. return(0);
  1340. }
  1341.  
  1342.